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Abstract 

A lifting body has been proposed as a candidate 
for the ACRV (assured crew return vehicle), which 
will serve as the crew rescue vehicle for Space Station 
Freedom. The focus of this work is on body surface 
definition, surface and volume grid definition, and 
the computation of inviscid flow fields about the vehi- 
cle at wind-tunnel conditions. Very good agreement 
is shown between the computed aerodynamic charac- 
teristics of the vehicle at Mach 10 and those measured 
in wind-tunnel tests at high Reynolds numbers. 

Introduction 

Currently within NASA, there is an ongoing ac- 
tivity to define candidate approaches for crew res- 
cue from Space Station Freedom. In-house NASA 
studies have shown that a station-based rescue ve- 
hicle will be required for assured crew return if the 
Space Shuttle is not operational, if there is an emer- 
gency that requires abandoning Freedom , or if there 
is a medical emergency that requires evacuation of 
the crew to Earth. A number of assured crew re- 
turn vehicle (ACRV) concepts that have been pro- 
posed for study over the next 2 years are shown 
in figure 1. The concepts range from blunt bod- 
ies, much like the Apollo shape, to a lifting-body 
shape derived from work initiated by NASA and the 
U.S. Air Force in the late sixties when studying the 
Dynasoar, HL-10, X-24, and M2-F2 vehicles. 

As part of the ACRV activity, the Space Systems 
Division at Langley Research Center has been study- 
ing the lifting-body concept. This proposed vehi- 
cle, which is currently envisioned to be about 25 ft 
long, will have the capability to return eight peo- 
ple to Earth from Space Station Freedom and land 
on conventional runways. The lifting-body concept 
is the most sophisticated of those being studied for 
the Space Station Freedom rescue mission. The com- 
plex aerodynamic shape of the vehicle can be seen 
in figure 1, along with the wings located outboard 
on the body and the small centerline vertical tail. 
Wind-tunnel tests indicate that the vehicle will have 
very good flight characteristics. The purpose of the 
work reported herein is to establish a computational 
fluid dynamics (CFD) capability to support system 
studies and experimental testing of the lifting-body 
concept. These capabilities would include a deter- 
mination of inviscid and viscous flow fields at a 
wind-tunnel condition and a comparison of computed 
aerodynamic characteristics of the vehicle with those 
determined through wind-tunnel testing. 

Symbols 

A inviscid flux Jacobian matrix 


A i , Bj, C i , Dj 

vectors of cubic-spline 
coefficients 


vectors of spline coefficients in 

C /,/b D /,fc 

streamline direction 

a, b 

inviscid flux volume weighting 
parameters 

a Xi ) b Xi , c Xi , d Xi , 
a 3/*> fy/i’D/i’ ^2/i 

cubic-spline coefficients 

C 

nondimensional speed of sound 

C A 

axial-force coefficient 

Cm 

pitching-moment coefficient 

C N 

normal-force coefficient 


derivative at inflection point 

D 

summation of flux differences 

dXj 

dt 

vector of first spline derivative 

rf 2 X, 

dt 2 

vector of second spline 
derivative 

E 

nondimensional total energy 
per unit volume 

E' 

dimensional total energy per 
unit volume, kg/m-sec 2 

E,F, G 

flux vectors in £, r), and C 
directions 

e 

nondimensional internal energy 
per unit mass 

e' 

dimensional internal energy 
per unit mass, m 2 /sec 2 

H 

flux vector normal to cell face 

h 

Cartesian flux vector 

/ 

index for grid points in cir- 
cumferential direction 

i,j,k 

indices of cell centers in £, rj. 
and C directions, respectively, 
or index of measured points 

K 

number of grid points in 
streamwise direction 

L 

normalizing reference length, 


m 

L/D 

lift-drag ratio 

M 

number of grid points in 
streamwise direction 


Moo 

free-stream Mach number 

V, W 

nondimensional velocities in 

M 

right inviscid Jacobian eigen- 


direction tangent to cell face 


vector matrix 

X 

vector of coordinate functions 

Mr 1 

left inviscid Jacobian eigenvec- 


vector of measured coordinates 


tor matrix 

X"(t) 

vector of second derivatives 

m 

number of measured sections 
in stream-wise direction 


nondimensional Cartesian 
coordinates 

N 

number of grid points around 
vehicle 

x(t),y{t) 

coordinates as a function of 
parametric variable 

n 

current time level 

x,y,z 

interpolated coordinates for 


number of measured points on 


streamwise smoothing 


kth section 

x', y', z' 

dimensional Cartesian coordi- 

n 

unit vector normal to cell face 


nates, m 

P 

nondimensional pressure 

Ui 

zth cross-sectional coordinates 

P' 

dimensional pressure, N/m 2 

a 

angle of attack, deg 

P 

dummy index for i, j, or k 

ri,r 2 

r 3 

grid concentration parameters 
concentration control parame- 

q 

dependent variable vector 

ter for streamwise direction 

R 

Reynolds number 

7 

ratio of specific heats 

S' 

velocity, m/sec 

A() 

spatial change in quantity, 

Si 

set of spline coefficients in 


Or - ()l 


streamwise direction 

<50 

time change in quantity, 

Sk 

set of all spline coefficients for 


( ) n+1 - ( r 


A;th section 


vector of allowable 

s 

streamwise parametric variable 


deviations 

t 

time, sec, or parametric 

e 

positive constant 


independent variable 

c 

computational coordinate in 

t 

intermediate parametric 


streamwise direction 


variable 

V 

computational coordinate 

te 

parametric variable value at 
inflection point 

X 

inviscid Jacobian eigenvalue 
matrix 

U 

ith value of parametric 
variable 


computational coordinate 
around configuration 

tn 

last value of parametric 
variable 

Se 

computational coordinate 
inflection value 

u 

nondimensional velocity in 

P 

nondimensional density 


direction normal to cell face 

P' 

dimensional density, kg/m 3 

u, v, w 

nondimensional Cartesian 
velocity components in x, y, 

a 

nondimensional cell-face area 


and 2 directions 

n 

nondimensional cell volume 

u , v 1 , w' 

dimensional Cartesian velocity 

Subscripts: 



components in x, y , and 2 
directions, m/sec 

L 

quantity at cell center to left 
of cell face 


2 


R 

oo 

Superscripts: 

n 

T 

+ 


* 

/ 


quantity at cell center to right 
of cell face 

free-stream value 


time level 
transpose 

quantity associated with 
positive eigenvalue 

quantity associated with 
negative eigenvalue 

latest available data 

dimensional quantities 


Geometry 

All computational fluid dynamics simulations 
require the generation of a grid that covers the flow- 
field domain and its boundaries. To obtain an accu- 
rate solution, the grid must be concentrated in re- 
gions with large gradients, and the grid must possess 
favorable mathematical characteristics, such as ac- 
ceptable grid skewness. 

For the flow about an aerospace vehicle such as 
the ACRV, a grid on the configuration boundary 
surface must be created and integrated into an overall 
grid. The vehicle surface grid is extracted from a 
numerical model obtained from either a computer- 
aided design (CAD) system or from measurement of 
a physical model. 

The numerical model of the ACRV vehicle was ob- 
tained from the measurement of a 5. 9-in- long wind- 
tunnel model. The measuring instrument obtains 
two coordinates around the model at fixed stations 
in the third coordinate. When the vehicle geometry 
is defined in a right-hand Cartesian coordinate sys- 
tem (fig. 2), the origin of the coordinate system is 
at the most forward point of the vehicle. The nega- 
tive 2 -axis extends down the vehicle, the y-axis is in 
the spanwise direction, and the .r-axis is in the verti- 
cal direction. The measured data base consists of x 
and y surface coordinates at 20 2 -coordinate stations 
(cross sections) down the vehicle. (See fig. 3.) 


or surface roughness. These irregularities are illus- 
trated in figure 4, which shows a plot of the original 
measured surface points at the 15th station along 
the body. The solid line represents a spline fit of 
the measured data points. Inspection of this plot 
does not reveal any noticeable surface imperfections. 
However, the expanded views of three surface sec- 
tions enclosed by the boxes labeled 1, 2, and 3 and 
shown in figures 4(b), 4(c), and 4(d), respectively, 
show surface irregularities present in the measured 
data. If the irregularities are allowed in the compu- 
tational grid, they will also appear in the numerical 
solution. Alternatively, the numerical model can be 
initially idealized by smoothing the measured surface 
coordinates. There may also be surface characteris- 
tics that are not to be included in a CFD solution. 
For instance, the vertical stabilizer is not included in 
the initial computations. 

Surface Smoothing 

To remove imperfections in the measured data, a 
cubic-spline smoothing algorithm (ref. 1) is applied 
to the coordinates at each 2 -coordinate cross section. 
The application of the algorithm starts with the cre- 
ation of a parametric independent variable t, which 
is the approximated arc length along a cross section. 
The x and y coordinates are considered to be func- 
tions x(t) and y{t). The functions are determined 
from the sets of parametric data 

{x uti }^ and {y u t t }‘i? 

where x t and y t are ith measured coordinates on a 
cross section, t, is the approximate arc length at the 
ith point, and n is the number of measured points on 
a cross section; n varies from cross section to cross 
section. The approximate arc length is defined by 

T 2 9 "I t/2 

^i+1 = ~ x i) T (Vi+1 ~ Vi) J T 

where t\ = 0. This computation starts on the bottom 
of the vehicle at the symmetry plane and continues 
around to the top of the vehicle at the symmetry 
plane. After the maximum approximate arc length 
is computed, the variable t is redefined to be the 
normalized approximate arc length by 


Measured coordinates have associated measure- 
ment errors, and there is the possibility of roughness 
on the model itself, which occurs in the manufac- 
turing process. An inspection of the plotted cross 
sections of the ACRV model revealed, from a com- 
putational standpoint, a number of serious irregular- 
ities that could be the result of measurement error 


ti = 


U_ 

tn 


The functions x(t ) and y(t) are represented in vector 
form 


X(t) = 


x( D 

vb) . 


3 



and the cubic-spline representation is 

X(i) = [a* + B i(t- ti) + c i(t- ti) 2 

-\ i=n — 1 

i=l 


where 


+ Bj(t — ti 


Aj = 


a Xi 

a Vi 


B i = 


L hi J 


( 1 ) 


C, = 


Cxj 

C Vi J 


D < = 


dx 

d 


Vi J 


The spline conditions are 


/ dX, 


V dt / ti+1 
/d 2 xA 

V " 2 L 


dX 


»+i \ 


dt J t 


i+l 


'fi 2 X 


8 + 1 | 

dt j 


+ 1 1 


where 


i = 1, 2, - ■ -, n — 2 
and where 0 < t <t n . The coefficients 

{Aj, Bj, Cj, 

are undetermined parameters whose values define 
X(t) in equation (1). The objective is to find the 
coefficients that minimize the integral of the second 
derivatives squared as follows: 


f 


[X"(t)] 2 dt 


subject to the constraint 


[X(tj) — X m (f i )] 2 < ^ 


8=1 


where 


fiX,- 


X m (tj) — 


(2) 


Ls/i. 


surface coordinates X m . The restated objective is 
to find the smoothest cubic spline that passes within 
the following bounds: 


X m (tj) — fiXj < XmCti) < Xi(t m ) + dXj 


where fiXj is the maximum deviation of the spline 
functions from the measured coordinates. The 
method of Lagrange multipliers from the calculus of 
variations is used to find the parameters 

Sfc = {Aj, Bi, Cj, 

where k is the fcth cross section in the streamwise 
direction. The solution algorithm for spline smooth- 
ing can be found in reference 1. The technique ex- 
ists in subroutine form (SUBROUTINE CSDSES) 
in the Langley math library. (See ref. 2.) A main 
program is written to read the measured coordi- 
nates, compute the approximate arc length, and call 
the spline-smoothing subroutine for each coordinate. 
It is necessary to provide the allowable deviations, 
{fiXj)A" fc , and the one constant e, which is an esti- 
mate of the sum of the normalized deviations squared 
in equation (2), is set equal to — 1. Figures 5 
and 6 show the application of the smoothing proce- 
dure to two cross sections from the ACRV vehicle 
(one near the front of the vehicle and one near the 
rear of the vehicle). Figures 5(a) and (b) and fig- 
ures 6(a) and (b) show the y versus t and x versus 
t smoothing; figures 5(c) and 6(c) show the x ver- 
sus y cross section. The algorithm is applied to each 
measured cross section. 

Surface Grid 

The technique described previously smoothes the 
measured data at each cross section in only one 
direction. The next step is to distribute points 
along the smoothed cross sections in a fashion that 
is desired for the final grid points. These points are 
then smoothed and interpolated in the streamwise 
direction for the surface grid points. The initial 
desired circumferential distribution of points on the 
ACRV vehicle is uniform around the vehicle near 
the front and concentrated at the wingtip as the 
wing emerges from the fuselage downstream. The 
streamwise distribution should be concentrated near 
the front of the vehicle. 


fiXi = 


8xi 

hi 


and E is a positive constant that specifies the extent 
of smoothing. The vector fiX is the allowable devi- 
ation of the spline function X(t) and the measured 


If N is the number of grid points around the 
body, and the index I denotes the ith grid point, 
a computational coordinate £ around the body is 
defined by 




(I- 1) 

(IV -1) 
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and £ is mapped into the normalized approximate 
arc length t with the variable t. Interpolated x and 
y coordinates are obtained from and equation (1) 
given t. In the streamline direction and up to the 
wing root, the desired grid spacing is uniform and 
t = £. At the wing root, the distribution becomes a 
bi-exponential that is defined by 

exp(riQ - 1 
e exp(r!) - 1 

(o < £ < £ e ; 0 < t < t e ) 


approximate arc length is computed for each value 
of I, s is normalized with respect to the maximum 
value. The same procedure that was used to repre- 
sent a cubic-spline approximation for X(t) is used to 
compute X(7, s) where 


X(I,s) 


x(I,s)' 

yl^s) 

z(/,s)_ 


and the cubic-spline representation is 


and 


X(I,s) 


A /,fc + B /,fc( s “ s I,k ) 


t — t e + (1 — t e ) 


exp (r 2 frrf) -1 
exp(T 2 ) - 1 


0£e < £ < 1; te < t < 1) 


The symbol £ e is the desired percentage of grid points 
around the bottom of the vehicle up to the wing lead- 
ing edge, and t e is the normalized approximate arc 
length along the same curve up to the wing leading 
edge. The constant Ti governs the amount of concen- 
tration at the wingtip, and T 2 is computed by using 
a Newton-Raphson iteration to assure continuity of 
(Dt/D£)!; e . The parameters £ e , te, and Ti vary in 
the streamwise direction from the wing root to the 
end of the wing. The computation of the surface grid 
is based on smoothing the previously defined data in 
the streamwise direction. The parametric form of the 
smoothed cross-sectional data is 

r - ,I=N,k=m 

\ x I,k i s I,k } /= i,fc=i 


r - i /= JV , fe=m 

{yi,k,si,ki I= i )jfc= i 


{ z k, s i,k}k=T 


where z & is the negative 2 coordinate of the /cth cross 
section, m is the number of cross sections, and s is the 
normalized approximate arc length in the streamwise 
direction defined by 


s I,k + 1 


( x I,k + 1 - x I,k ) 2 + {yi,k+l ~ Vl,k ) 2 


+ ( z k + 1 — z k ) 2 


1/2 


+ s /,fc 


S I, l = 0 


+ C I,k( s - s I,k ) 2 


+ D /,fc( s “ s I,kf 


1 1=N,k=m 


I=l.k=l 


( 3 ) 


where 


ax I,k 

_ 

^ x i,k 

a yi,k 

II 

byi,k 

. az I,k . 


. b z i,k . 

Cx I,k 


d*i,k 

a Vl.k 

II 

IQ 

dyi,k 

. az I,k . 


. b zi,k - 


The procedure described above is used to compute 


Si 


A I,k > ^ I,k , 


I=N,k=M 
/= l,fc=l 


where 


<5X i tk = 


bxi,k 

6 yi,k 


If M is the total number of grid points along the 
body in the streamwise direction, and if the index k 
is the fcth grid point, a computational coordinate ( 
along the vehicle is defined by 


c = 


fe - 1 
M - 1 


and £ is mapped into the normalized approximate 
arc length s. A concentration of grid points near 
the front of the vehicle is achieved by using the 
transformation 


The computation of s /j. starts at the nose of the ve- 
hicle and ends at k = m. Again, after the maximum 


exp(T3<f) — 1 
exp(T 3 ) - 1 
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The constant 1^ determines the amount of concen- 
tration. A surface grid 

[x(7, K),y(I, K), z(I, 

is obtained by evaluating equation (3) and with a 
given value of s. 

Starting with the measured data, a surface grid 
for the ACRV vehicle where N = 81 and M = 125 is 
shown in figure 7. This boundary surface grid is then 
used in the computation of the volume grid about the 
vehicle on which Euler flow is simulated. 

Volume Grid 

The initial flow-field computations about the 
ACRV are obtained at a high Mach number and a 
high angle of attack. Under these conditions, the 
bow shock lies close to the body on the windward 
side and far from the body on the leeward side. The 
need for the volume grid to extend beyond the outer 
limits of this shock dictates the structure of the vol- 
ume grid. 

The grid topology chosen for the ACRV flow field 
is a dual-block topology similar to that described in 
reference 3. In this topology, there is an inner block 
that begins from a singularity line just upstream 
of the wing root. A second block surrounds the 
forward part of the vehicle and the inner block. 
The general representation of the two-block topology 
is shown in figure 8. In reference 3, the lifting 
surfaces on the configurations have sharp leading 
edges, whereas the ACRV wing has a rounded leading 
edge. Consequently, the outer-boundary surface of 
the inner-grid block is displaced forward of the wing 
leading edge so that a C-type grid is formed spanwise 
about the wing. 

The advantage of a two-block topology over a 
one-block topology is that the flow-field region can 
be adequately covered with fewer points. That is, 
to have enough points to cover the region between 
the wing and symmetry plane, a single block would 
require many more points in the forward region than 
is necessary. The disadvantages are the additional 
program complexity in the solution code and the 
numerical complications at the singularity line. 

Inner-Block Volume Grid 

The inner block is bounded by the vehicle surface, 
the outer-boundary surface of the inner block, the 
bottom-symmetry surface, the top-symmetry sur- 
face, the singularity line, and the downstream sur- 
face (fig. 9). The outer-boundary surface is defined 


by creating an analytical curve, translating the curve 
to the singularity line, and rotating it around the 
vehicle. It is then stretched in the spanwise direc- 
tion and skewed upward as it passess over the top 
of the vehicle. This process is similar to the sur- 
face boundary generation that is described in refer- 
ence 3, with the addition of the stretching and skew- 
ing. The interior grid and the distribution of grid 
points on the outer-boundary surface, the symmetry- 
plane surfaces, and the downstream surface are com- 
puted by using a three-dimensional version of the 
two-boundary grid generation (TBGG) technique de- 
scribed in reference 4. The grid points are concen- 
trated at the vehicle surface, and near orthogonal- 
ity is maintained. Figure 10 shows the inner grid 
on the downstream boundary surface. As the line 
of singularity is approached from downstream, the 
distribution becomes more uniform and matches the 
distribution upstream of the singularity line. 

Outer-Block Volume Grid 

The grid on the outer block surrounds the forward 
part of the fuselage and the grid on the inner block. 
The six surfaces that form the outer block are (1) the 
forward fuselage and outer-boundary surface on the 
inner block, (2) a far-held boundary surface, (3) the 
symmetry plane on the bottom side of the vehicle, 

(4) the symmetry plane on the top side of the vehicle, 

(5) a polar singularity line from the (0,0,0) point of 
the vehicle to the far-held boundary surface, and 

(6) a downstream surface that extends the inner 
block to the far-held boundary. 

The hrst computation for the outer block is the 
far-held boundary surface, which is similar to com- 
puting the outer-boundary surface for the inner-grid 
block. The characteristics of the far-held boundary 
are that it is relatively close to the vehicle on the bot- 
tom side and relatively far away on the top side to 
be just outside of the bow shock. Also, it should ex- 
pand away from the wings in the spanwise direction 
to capture the shock phenomenon. 

Instead of using a single analytical curve and ro- 
tating it around the vehicle, two curves are blended 
in the streamwise direction to create another curve 
that is translated to a point about four shock standoff 
distances from the nose and rotated about the vehi- 
cle while being stretched in the spanwise and vertical 
directions. The two blended curves are chosen such 
that one is suitable forward and the other is suit- 
able rearward. Grid points are distributed rearward 
and around the outer boundary to be approximately 
across from the corresponding points on the surface 
of the opposite side. A side view of the far-field sur- 
face is shown in figure 11. 
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The general approach for computing the outer- 
block grid is transfinite interpolation (ref. 5), given 
grid point information at block boundaries and in- 
termediate surfaces. The outer-block computation is 
divided into two parts: (1) the forward part from 
the nose singularity line back to a surface extend- 
ing from the singularity line around the fuselage to 
the far-field boundary (singularity-line surface), and 
(2) the rearward part extending from the singularity- 
line surface to the downstream surface. 

The forward part of the outer block is obtained 
by first computing grid points on the bottom sym- 
metry plane, on the top symmetry plane, and on 
a spanwise surface nearly orthogonal to the sym- 
metry plane. The TBGG algorithm is applied to 
concentrate grid points near the vehicle surface and 
maintain near orthogonality. Figure 12 shows the 
three surfaces. Transfinite interpolation' with linear- 
exponential blending functions is used to fill in the 


interior. The rearward part on the outer block is 
computed by using transfinite interpolation with La- 
grangian blending functions. (See ref. 6.) This 
method allows first- derivative continuity in grid spac- 
ing and grid-curve direction except in the vicinity of 
the singularity line. A concentration of grid points in 
this region is used to overcome this deficiency. Fig- 
ure 13 shows grid points in the symmetry plane and 
two spanwise surfaces in the streamwise direction. 

Solution Technique 

The solution technique is based on the Langley 
Aerothermodynamic Upwind Relaxation Algorithm 
(LAURA) code discussed in reference 7; LAURA is 
a viscous code that has been modified to compute 
inviscid flow fields. A brief description of the LAURA 
code, which is an upwind finite-volume algorithm, 
and the implementation of the inviscid boundary 
condition are given in this section. 


Governing Equations 

The integral form of the governing equations can be expressed as 



q t dQ + 



h ■ n dcr = 0 


( 4 ) 


where q t is the time rate of change of the dependent variable vector at a cell center, h is the Cartesian flux 
vector, and n is the outward unit normal vector to a cell face. The cell volume and cell-face area are represented 
by U and a, respectively. The quantity h • n is the flux vector normal to a cell face. The vector h can be 
written in terms of its components in the x , y, and 2 directions as 


h = E x + Fy + G z (5) 

where 

Q =[p,Pu,pv,pw,E] r (6) 

E = \pu, pu 2 + P, puv, puw , (E + P)u] T (7) 

F = [pv, puv , pv 2 + P, pvw, ( E + P)v] T (8) 

G = \pw, puw, pvw, pw 2 + P,(E + P)w) T (9) 


All quantities have been nondimensionalized as shown in the following equations: 


P = 
P = 


x = 


d U S' 

Poo °oo 


V = 


S'oo 


w 


IV 

SL 


p' 


e = 


(S'o 


E = 


E' 

P'^S’oc) 2 


V y= u z= L> 


The value L' is a user-defined reference length. The equation set is closed by using the idea! gas law. 
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Expressing equation (4) in finite-volume form for a single six-sided cell in the computational domain gives 





a i-l/2 ,j,k - H *+l/2 ,j,k a i+l/2,j,k 


+ H i,j-l/2,fc <r ij'-l/2,fc H i,j+l/2,fc a i,j+l/2,k + H iJ,fc-l/2 <7 ij,fe+l/2 


1/2 flr i,j,fc+l/2j 

A shorthand index notation enables equation (10) to be expressed as 


c l' St 

= \ o 


o, „ E [(HV ) p _ 1/2 - (HV) p+1/2 ] 

hJ,k p=ij,k 


(10) 


( 11 ) 


where <5q is the change in q per time step, H is the inviscid flux normal to a cell face, and 6t is the time step. 
The lowercase subscripts in equations (10) and (11) represent cell-centered values, unless offset by one half, in 
which case they represent values at the center of a cell face. A diagram of the indexing is shown in figure 14. 
The equations used to compute the cell volume, the cell-face area, the time step, and the metrics are found in 
reference 7. 

Because the flux normal to a cell face is a function of information on either side of the face (left or right), 
H in equations (10) and (11) is expressed as 

H*(qL,q*) - i[aH*(q L ) + AH*(q fl )] - ^D*(q L ,q*) (12) 

where L and R denote the indices of the cell centers to the left and right, respectively, of the face being 
evaluated. The functions a and b are flux weighting parameters defined in terms of cell volumes as 

2 2 fii 

Qr + &L ^r + Q,l 

The parameters a and b lessen the effects of grid stretching, and their formulation is empirical. (See ref. 8.) 
The variable D* in equation (12) is a summation of flux differences. The value D* can be determined from 
the solution of a Riemann initial-value problem. Several solution procedures for D* have been developed. The 
approach adopted in this work is attributed to Roe (ref. 9). 

The first-order upwind-biased flux is computed by using the Roe flux-difference splitting (ref. 9). For the 
first-order flux, the flux difference across the cell face comes from the solution of a Riemann initial- value problem 
posed by Roe. The initial data for the problem are found at the two cell centers to which the face is common. 
The flux difference is split into two parts, one associated with the positive and one with the negative eigenvalues 
of the unsteady inviscid Jacobian matrix. The speed and direction of information propagating toward a cell face 
are directly proportional to the absolute magnitude of the eigenvalues and sign of the eigenvalues, respectively. 
That is, negative eigenvalues send information from the right toward a cell face, whereas positive eigenvalues 
send information from the left; how quickly this information travels is related to the absolute magnitude of the 
eigenvalue. Roe’s method correctly interprets this wave motion relative to a cell face and computes an inviscid 
flux that correctly models the physics of the flow. The inviscid flux for a supersonic flow with all positive 
eigenvalues is made up solely of information from the left, and a supersonic flow with all negative eigenvalues 
is made up solely of information from the right. If the flow has eigenvalues with mixed signs, the inviscid flux 
is comprised of information from the right and left. 

Roe’s method gives D as the exact solution to the following approximate Riemann problem: 


qt + Aq*. = 0 


(13) 
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The change in flux across a cell face for this problem can be expressed as 


D (qL,qj?) = AH 7 = - H/(q L ) = MAM 2 (q L , <Ir){<IR - qi) ( 14 ) 


where the flux vector Jacobian A is expressed in terms of its right (M) and left (M _1 ) eigenvectors and a 
diagonal matrix of its eigenvalues (A). The matrices A, A, M, and M -1 are given in reference 7. Representing 
the positive and negative contributions of the eigenvalues as shown below 


- A +1 A I X- A ~I A I 

2 2 


(15) 


enables D(qx,q#) to be written as 


D(q L , q*) = (MA + M~ 1 - MA'M-^qfl - q L ) 


(16) 


or 


D*(q L , q*) = M|A|M -1 Aq* = |A|Aq* (17) 

When equation (16) is substituted into equation (12), the first-order inviscid flux normal to a cell face is 

H}(qL,q*) = \ [oHj(qL) +6H}(qfl) ~ |A|(q fl) q L )Aq*] (18) 

The inviscid flux shown in equation (18) can be thought of as being composed of a second-order approximation 
to the inviscid flux at a cell face (first two terms) minus a dissipation term (remaining term). If this dissipation 
is not included, the algorithm is equivalent to a centrally differenced algorithm. 

For second-order-accurate solutions obtained by using Yee’s STVD approach (ref. 10), D* in equation (12) 
can be shown as 


D *(qi>qi?) = M 2 |A| 2 [M 2 x Aq^ - minmod (M 1 1 Aqi,M 2 1 Aq 2 ,M 3 1 Aq 3 )] (19) 


The subscript 2 references the face at which the flux is being computed — 1, the face behind, and 3, the face 
ahead. The minmod function compares differences in characteristic variables at these locations and chooses 
the smallest in absolute magnitude if the signs of the values are the same, or zero if the signs are different. 
The first term in equation (19) is the first-order term from equation (17), and the second is a correction that 
makes D* second-order accurate if the signs of the differences in characteristic variables being evaluated in the 
minmod function are the same. 


If variables with the superscript * contain information referenced to the i,j, k cell center, they are linearized 
by the following equation: 


K* = K n + 



*5q i,j,k 

i,j,k 


The variable K is a dummy variable which represents 


F V» 


O xro hi 41 


Vifii n rr 


9 


Upon substituting equation (12) into equation (11) and performing the requisite linearizations, the governing 
equations can be written as 


7 + 


2 n) ijk ^ [(I a I ct Wi/2 + (|a|<t)] 1/2 }**j,fc 

l ’ 3 ' k P=i,j,k 3 

]T (UaH?( qL )+bH?( qjt )-D n ( qL ,q R )}}a p _ 1/2 

t,3,k P=i,j,k 

{[aH?(q L ) + 6H?(q*) - D"(q L , q*)]} a P+ 1 / 2 ) 


( 20 ) 


Boundary Conditions 

This algorithm requires that dependent variables 
be specified both on the surface and at “pseudo” cells 
located one-half cell below the vehicle surface. In 
figure 14, (a) denotes the plane of cell-centered values 
one-half cell above the body surface, (b) denotes 
the plane of cell-centered values one-half cell below 
the body surface, and the X denotes the location 
of points on the surface. The values at (b) are 
required for the computation of the first- and second- 
order fluxes associated with the methods of Roe 
and Yee, respectively. The values used at (b) were 
based on the values at (a). The Cartesian velocity 
components at (b) were computed such that the 
following conditions were satisfied: 

U 6 = -U a 
v 6 = v a 
w 6 = w a 

where U is the contravariant velocity in a direction 
normal to the body surface, and V and W are 
contravariant velocities in a plane tangent to the 
body surface. First-order extrapolation was used to 
compute the internal energy at (b), and the density 
at (b) was computed by using 


'(7 - 1 )pe 


'(7 - 1 )pe 

L P 7 J 

a 

L P 7 


Values on the wall were determined such that 
surface tangency is enforced. (See ref. 11 .) In the 
method outlined in reference 11 , extrapolation is 
used to determine initial wall values, and a wave 
correction is then performed on these extrapolated 
quantities to determine final wall values and to sat- 
isfy tangency. In general, just extrapolating values to 
the surface will not satisfy tangency. The wave cor- 
rection is shown below, where c represents the wave 


corrected values and a represents the extrapolated 
values: 

Pa_Pc 

7 7 

Pa Pc 


u c = o 

V a - V c 
W a = W c 

At the remaining boundaries, pseudo cells were also 
employed. For the outflow boundary, a supersonic 
outflow boundary condition was used; along the axis 
singularity and symmetry planes, a reflective bound- 
ary condition is used. 

Flow-Field Computations and Analysis 

Preliminary trajectory analysis indicates that the 
ACRV vehicle will fly at approximately 25° angle 
of attack during the hypersonic portion of its flight. 
Also, a considerable amount of aerodynamic testing 
will be done in the Langley 31-Inch Mach 10 Tunnel. 
(See ref. 12.) Therefore, Moo = 10 and a = 25° 
have been chosen as the nominal conditions for the 
following flow-field solution. Also, solutions have 
been obtained at angles of attack of 15°, 20°, 30°, and 
35° to cover the range of available experimental data. 
Previous analyses (refs. 13 and 14) have shown that 
a perfect-gas analysis is appropriate when comparing 
computations with data from this tunnel. Therefore, 
the ratio of specific heats, 7 , was taken to be 1.4 for 
these computations. 

These computations were made on the CRAY-2 
computer at Ames Research Center’s Numerical 
Aerodynamic Simulation (NAS) system and on the 
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CRAY-2 computer at Langley Research Center with 
the body and grid defined in this paper. This grid of 
584 496 points translates into a code requirement of 
approximately 40 megawords of memory. A complete 
flow-field solution required about 10 hr of CPU time 
on the NAS CRAY-2 and about 8 hr of CPU time 
on the Langley CRAY-2, which is approximately 1.2 
times faster than the NAS machine. 

Figure 15 shows a pressure contour plot in the 
symmetry plane for the nominal condition. The ex- 
tent of the bow shock and the location of the canopy 
shock can be seen. The apparent smearing of the 
bow shock approximately half way down the body, 
and most evident in the windward-symmetry plane, 
is the result of the shock leaving the “outer” grid, 
which is closely aligned with the shock, and enter- 
ing the “inner” grid, which is highly skewed rela- 
tive to the shock. Nondimensional surface pressures 
in the symmetry plane are shown in figure 16. The 
windward-side pressure distribution is typical for this 
shuttle-like surface. On the leeward side, the flow 
expands rapidly around the nose, compresses when 
it encounters the canopy, and then overexpands and 
recompresses on the flat upper surface. These plots 
indicate that the smearing of shock in the inner grid 
has no effect on surface pressure. Additional repre- 
sentations of the gross characteristics of the external 
flow field are shown in figures 17 and 18. Figure 17 
shows pressure contour plots in three planes along the 
vehicle. The pressures lie in a computational plane 
normal to the axis of the vehicle. However, when 
the computational plane is transformed into Carte- 
sian space, the planes are not flat but are curved 
three-dimensional surfaces as they appear in this fig- 
ure. Here, black represents the lowest pressure and 
white the highest. The extent of the bow shock is 
clearly visible along with the rapid expansions about 
the chines along the side of the vehicle. For each 
plane shown, the complete computational plane is 
not displayed, and the visual quality of the picture 
is therefore improved. In figure 18, Mach number 
contour plots are shown in planes that correspond to 
those used in figure 15. The Mach number contours 
reveal more detail of the flow in the leeward region. 
In the second plane, the beginning of a cross-flow 
shock can be seen near the upper symmetry plane 
near the wall. This shock can be seen quite clearly 
in the third or outflow plane. Also, there appears to 
be a small wing-body shock interaction in the third 
plane, as indicated by the inflection point in the bow 
shock near the wingtip and by the high pressure at 
the wingtip in the third plane of figure 17. 

Figures 19 to 21 are all similar in that they show 
meridional distributions of nondimensional pressures 


along with the corresponding surface contour at three 
z locations along the body. The pressure is plotted 
as a function of y, so that the location of the pressure 
on the surface can be noted. The pressure is plotted 
on a log scale, so that pressure fluctuations on the 
leeward side can be seen. As expected, there is a 
slight rise in pressure when going from the centerline 
outboard at all locations. On the leeward side, 
pressure behaves as expected — the flow expands and 
compresses with corresponding changes in surface 
slopes. In each figure, the symbols on the surface 
contour plot represent the distribution of grid points 
on the surface. 

Figure 22(a) shows an oblique view of the highly 
three-dimensional exit surface of the volume grid. 
The corresponding pressure contour plot is shown in 
figure 22(b). 

Experimental data were taken in the Langley 
31-Inch Mach 10 Tunnel at angles of attack be- 
tween 15° and 30° at Reynolds numbers per foot of 
0.5 x 10 6 and 2 x 10 6 . The location of the moment 
center and the convention for positive aerodynamic 
coefficient are shown in figure 2. A comparison of the 
computed and measured axial-force coefficient C 4 is 
shown in figure 23. For a body such as the ACRV, the 
predominant effect of viscous forces should appear in 
Cj^. As shown in figure 23, the values of C 4 approach 
the computed inviscid limit as Reynolds number is 
increased. However, for the normal-force coefficient 
Cyv, shown in figure 24, there are apparently little if 
any viscous effects. 

A comparison of the computed and measured 
pitching-moment coefficient C m is presented in fig- 
ure 25. The measured data indicate little if any ef- 
fect of Reynolds number on C m , and the computed, 
inviscid values of C m are generally in good agree- 
ment with the measured data. The trim angle of at- 
tack inferred from both sets of data would be almost 
identical. There is a significant impact of Reynolds 
number on the lift-drag ratio, L/D, as indicated in 
figure 26. However, both the computed L/D and 
measured L/D show that, regardless of viscous ef- 
fects, the maximum L/D occurs at the trim angle of 
attack. 

Conclusions 

In this paper, the tools necessary to support sys- 
tem studies and hypersonic aerodynamic testing of a 
lifting-body concept have been established. It has 
been demonstrated that the Langley Aerothermo- 
dynamic Upwind Relaxation Algorithm (LAURA) 
code, modified for inviscid comput at 10ns, has 
versatility and robustness to handle the flow about 


11 



complex geometries at hypersonic speeds and high 
angles of attack. 

A methodology has been established for con- 
structing a surface definition that meets the require- 
ments for smoothness and continuity required for a 
quality computational fluid dynamics analysis. Also, 
a multiblock approach combined with an algebraic 
grid generation technique has been shown to produce 
a useful volume grid about a complex geometry. 

The aerodynamic coefficients generated by using 
the grid and flow solver described herein are in 
very good agreement with the measured data at 
high Reynolds number conditions from wind-tunnel 
tests. Although overall vehicle characteristics, such 
as aerodynamic coefficients, are poor indicators of 
the performance of a flow solver or a grid, they do 
indicate that the general characteristics of the flow 
field are being properly modeled. A more detailed 
analysis of the capability of the flow solver and grid 
will have to await the availability of detailed surface 
data. 

NASA Langley Research Center 
Hampton, VA 23665-5225 
June 24, 1991 
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(a) Measured surface definition. 






















Figure 13. Composite of grids in symmetry plane and two spanwise surfaces. 
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Figure 23. Computed and measured axial-force coefficients. = 
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Figure 24. Computed and measured normal-force coefficients. M 0 0 = 10 and 7 
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Figure 26. Comparison of computed and measured lift-drag ratios. 
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